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We show that the van der Waals-like phase separation instability of a driven granular gas at 
zero gravity, previously investigated in two-dimensional settings, persists in three dimensions. We 
consider a monodisperse granular gas driven by a thermal wall of a three-dimensional rectangular 
container at zero gravity. The basic steady state of this system, as described by granular hydrody- 
namic equations, involves a denser and colder layer of granulate located at the wall opposite to the 
driving wall. When the inelastic energy loss is sufficiently high, the driven granular gas exhibits, 
in some range of average densities, negative compressibility in the directions parallel to the driving 
wall. When the lateral dimensions of the container are sufficiently large, the negative compressibility 
causes spontaneous symmetry breaking of the basic steady state and a phase separation instability. 
Event-driven molecular dynamics simulations confirm and complement our theoretical predictions. 

PACS numbers: 45.70.Qj 



I. INTRODUCTION 

Rapid flows of granular materials are widespread in 
nature and technology. Losing kinetic energy to micro- 
scopic degrees of freedom of the grains in grain collisions, 
the granular flows are intrinsically far from thermal equi- 
librium and therefore exhibit a host of pattern forming 
instabilities [J Quantitative modeling of granular 
flows remains challenging, and pattern forming instabil- 
ities can help by providing sharp tests to the models. 
In this work we focus on one pattern forming instability 
that develops in granular gas: an assembly of inelasti- 
cally colliding hard spheres. The only dissipative effect 
in the particle collisions that we will take into account 
is a reduction in the relative normal velocity of collid- 
ing particles, accounted for by a (constant) coefficient of 
normal restitution r. We will assume nearly elastic col- 
lisions, 1 — r <C 1, and small or moderate gas densities. 
As shown in many previous studies [3|, U, [5|, [fj, [7|, [8| , these 
restrictions make it possible to use equations of granular 
hydrodynamics. 

The phase separation instability, that will be in the 
focus of our attention here, was originally predicted from 
hydrodynamic equations and then observed in molecu- 
lar dynamic (MD) simulations in a two-dimensional (2D) 
setting: a monodisperse gas of inelastically colliding hard 
disks at zero gravity, confined in a 2D rectangular box 
and driven by a side wall that vibrates with a high fre- 
quency and small amplitude d, EH EH EH EH EH EH • 
The basic steady state of the 2D system is the stripe state 
[l7l EH : a stripe of a denser and colder gas located at 
the wall opposite to the driving wall. At sufficiently high 
energy loss, and within a certain ("spinodal") interval 
of grain area fractions, the stripe state becomes unsta- 
ble with respect to small density perturbations in the 



lateral direction, unless the lateral container size is too 
small [H EH EH] ■ Within a broader binodal, or coexis- 
tence, interval, the stripe state is metastable [TH EH- In 
both cases one finally observes, usually after a coarsen- 
ing process, a granular "drop" coexisting with "vapor", 
or a granular "bubble" coexisting with "liquid" , along 
the wall opposite to the driving wall [TH EH EH ■ This 
remarkable far-from-equilibrium two-dimensional (2D) 
phase separation phenomenon is in many ways similar to 
the gas-liquid transition as described by the classical van 
der Waals equation of state, but the role of temperature 
is now played by the inelastic energy loss, see below. The 
basic properties of the phase separation in 2D are qualita- 
tively captured by an effective one-dimensional granular 
hydrodynamic model, suggested in Ref. [T(|. Recently, 
the studies of the phase separation in 2D have been ex- 
tended to an annular geometry [l9| . 

The present work predicts a similar phase separation 
phenomenon in three dimensions (3D). By extending the 
previous treatments to 3D we are breaking ground for a 
future investigation of this phase separation phenomenon 
in reduced gravity experiments. The paper is organized 
as follows. In Sec. |TT]we employ a hydrodynamic descrip- 
tion to describe the "layer state" (the basic state of the 
system), to compute the spinodal balloon and the binodal 
asymptote of the system, and to determine the critical 
lateral dimensions of the container for the phase separa- 
tion to occur. As this hydrodynamic description will be 
dealing only with steady states with a zero mean flow, the 
corresponding theory may be called hydrostatic. In Sec- 
tion IIIII we report a series of event-driven MD simulations 
of this driven granular system and compare the simula- 
tion results with the hydrostatic theory predictions. Sec- 
tion IIVI summarizes our results, discusses possible mor- 
phologies of phase-separated states and briefly mentions 
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some extensions of the model. 



II. GRANULAR HYDROSTATICS: THE LAYER 
STATE AND THE PHASE SEPARATION 

A. The density equation 

Let N hard spheres of diameter d and mass m = 1 
move, at zero gravity, inside a rectangular container with 
dimensions L x , L y and L z . The spheres collide inelasti- 
cally with a constant coefficient of normal restitution r. 
For simplicity, we neglect the rotational degree of free- 
dom of the particles. Let one of the container walls per- 
form high-frequency and small-amplitude vibrations in 
the x-direction. We assume that the vibration amplitude 
is much less than the mean free path of the particles at 
the driving wall, while the vibration frequency is suffi- 
ciently high. In this case one can treat the vibrating 
wall as effectively immobile and prescribe a constant gas 
temperature To at this wall 20]. For simplicity, particle 
collisions with all other walls of the container are consid- 
ered elastic. The energy transferred from the "thermal" 
wall to the granulate dissipates in the inter-particle col- 
lisions, and we assume that the system reaches a time- 
independent state with a zero mean flow. In the nearly 
elastic limit, 1 — r <C 1, and for small or moderate parti- 
cle densities, one can safely use granular hydrodynamic 
equations. For a zero-mean-flow steady state these are 
reducible to two hydrostatic relations: 



V ■ [kVT(i-)] 



p = const . 



(1) 



where the gas pressure p = p(n, T), heat conductivity k — 
«(n, T) and energy loss rate due to inelastic collisions I — 
I(n,T) depend on the particle number density n(r) and 
granular temperature T(r) of the gas. We will employ the 
equation of state of Carnahan and Starling [2l| and the 
widely used semi-empiric transport coefficients derived 
from kinetic theory in the spirit of Enskog approach [22| : 
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and v — Trd 3 n/6 is the local value of the solid fraction of 
the grains. Let us rescale all the coordinates by L x and 
introduce the rescaled inverse density u(r) = n c /n(r), 
where n c = \/2/d 3 is the crystal packing density in 3D. 
The rescaled coordinates x, y and z now run between 



while Eqs. can be transformed, after some algebra, 

into a single equation for the inverse density u(r): 



V • [F(u)Vu] = AQ(u) . 
where F(u) = A(u)B(u), 
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while A = (7r/3)(l — r 2 )(L x /d) 2 is the dimensionless 
hydrodynamic inelastic loss parameter. The boundary 
conditions for Eq. (JSJ) are the zero heat flux conditions 
d n u = at all walls except the thermal wall x = 1, and 
the condition d y u = d z u = at x = 1. The latter con- 
dition follows from the constancy of the temperature at 
the thermal wall [201 ] . combined with the constancy of 
the pressure in a steady state. As the total number of 
particles N is conserved, we obtain 
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where / = (n) /n c is the average volume fraction of the 
granulate, and (n) — N / (L x L y L z ) is the average number 
density of the particles in the container. The nonlinear 
partial differential equation (PDE) ([5]), together with the 
boundary conditions and the normalization condition (|7|). 
determine all possible steady state density profiles, gov- 
erned by two dimensionless parameters / and A. The 
density profiles are independent of To- Importantly, at 
large u the function Q(u) decreases with an increase of u. 
This implies that the steady state solution of Eq. (0) is 
non-unique [23| which paves the way to phase separation 
and coexistence, as in the 2D case. 



B. The layer state, spinodal balloon and binodal 
asymptote 

The simplest solution of Eq. © u — U (x) is laterally 
symmetric, that is independent of y and z. This is the 
"layer state" , and it is fully determined by the following 
equations: 



[F(U)U']' = AQ(U), U'(0) = 0, 



dx 



= f, (8) 



/o U(x) 

where the primes denote the cc-derivatives. Figure [T] de- 



zero and 1, A y = L y /L Xl and A z — L z /L Xl respectively, picts an example of the rescaled density profile of the 
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layer state obtained by solving Eqs. fSJ numerically The 
hydrostatic density profile agrees with a late-time profile 
observed in our MD simulations described below. 




FIG. 1: The rescaled density profiles obtained from the 
hydrostatic equations (line) and from MD simulations (cir- 
cles). The dimensionless parameters are A = 5 x 10 and 
/ = 0.02317, the MD simulation parameters are L x — 500d, 
L y = L z = 50d, N = 40960 and r = 0.9904. The inset 
shows a typical snapshot of the system at the steady state as 
observed in the MD simulation. 
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FIG. 2: Upper panel: the rescaled steady state granular pres- 
sure P versus the grain volume fraction / for A = 199.3 and 
A = 5 x 10 3 . Shown by the solid line is the spinodal balloon, 
inside which the effective lateral compressibility of the gas is 
negative. The borders /i and f% of the spinodal interval are 
determined from the condition (dP/df) A = 0. Lower panel: 
a zoom-in at the P(f) dependence for A = 5 x I0 3 . 

Having found the density profiles at different A and /, 
we can compute, with the help of Eq. @, the rescaled 
pressure of the layer state P = p/(n c T ). As the steady- 
state pressure is constant throughout the system, one can 



compute it at the thermal wall x = 1, where the temper- 
ature T — T is prescribed [2(J. We obtain, therefore, 



l + 4G[t/(l)] 
U(l) 



Two typical P(/) curves for different A are shown in 
Fig. H At small A (exemplified by A = 199.3) the bulk 
energy loss is not very important, and P(f) is monotone 
increasing with /. At sufficiently large A (exemplified by 
A = 5 x 10 3 ) there is an interval of the volume fractions 
/i(A) < / < /2(A) where the rescaled pressure P(/) de- 
creases with an increase of /. Therefore, the effective 
compressibility of the gas in the lateral directions is neg- 
ative there. The lower panel of Fig. [5] shows a blowup of 
the negative compressibility region at A = 5 x 10 3 . The 
borders of the negative compressibility region are deter- 
mined by the condition (dP/df) A = 0. By joining the 
spinodal points /1 and (separately) fa at different A, we 
obtain the spinodal balloon of the system in the (P. /) 
plane (Fig. [j), or in the (A,/) plane (Fig. [3]). As A de- 
creases, the spinodal interval fi (A) < / < /2(A) shrinks 
into a point, as in the 2D case [10, [3 • This is the critical 
point of the system (P c , f c ), or (A c , f c ). At A < A c P(f) 
is monotone increasing. 

A negative lateral compressibility implies that, within 
the spinodal balloon, the layer state [a ID solution of 
the steady state equation |(SJ)] is unstable with respect 
to small-amplitude long-wavelength perturbations in one 
or both lateral directions. Similarly to the well-studied 
2D setting [l(J EH there is also a binodal (or coexistence) 
line, originating from the intervals of area fractions where 
the layer state, although linearly stable, is nonlinearly 
unstable (that is, metastable). The two branches of the 
binodal line in the (A, /) plane merge at the same crit- 
ical point (A c ,/ C ). The asymptote of the binodal line 
in a close vicinity of the critical point, |/ — f c \ <C f c 
and < A — A c <C A c , can be readily established, cf. 
Refs. [r| EH- Indeed, in the close vicinity of the crit- 
ical point P(f) is describable, at fixed A, by a cubic 
parabola in / — f c without a quadratic term. As a re- 
sult, one can find, at fixed A, the two points /_ and /+, 
belonging to the binodal line, from the simple relations 
P(f-A) = P(f + , A) and /_ + /+ = 2/ c . The resulting 
binodal asymptote is depicted in Fig. 3. 

Unfortunately, this simple asymptote cannot be con- 
tinued beyond the close vicinity of the critical point. The 
form of the binodal line far from the critical point has not 
yet been derived from granular hydrodynamics, neither 
in 2D, nor in 3D. Such a derivation would require a non- 
perturbative solution of the nonlinear PDE ([5]). Most 
likely, this can only be done numerically. Note that the 
"Maxwell construction", suggested in Ref. [r| for the 
binodal line in 2D, is valid only in a close vicinity of the 
critical point, where it is reducible to the two simple re- 
lations P(/_,A) = P(/+,A) and /_ + f+= 2f c . The 
reader is advised to consult with Ref. [15| for a more 
detailed discussion of this issue. 
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FIG. 3: The spinodal balloon of the system (the solid line) 
and the asymptote of the binodal line in a close vicinity of 
the critical point (the dashed line) on the plane /, A -1 ' 2 (the 
upper panel) and on the plane /, P (the lower panel). The 
insets zoom in on a close vicinity of the critical point. 



C. The critical value of lateral dimensions: a 
marginal stability analysis 

When the dimensionless parameters A and / are within 
either spinodal, or binodal balloon, a steady state with 
a broken lateral symmetry should develop. However, the 
phase separation demands a sufficiently large lateral size 
of the system. It will be suppressed by the lateral heat 
conduction if the lateral aspect ratios A y and A z are 
both less than a critical value A C (A, /). By analogy with 
2D, see Refs. [I [H, [H, H], we can determine A C (A, /) 
from a marginal stability analysis. Indeed, let A y and A z 
be less than A C (A, /), so the layer state is linearly stable, 
because of the lateral heat conduction, even within the 
spinodal balloon. Increasing A y and/or A z slightly be- 
yond A C (A, /), one should observe a (weakly) phase sep- 
arated steady state that bifurcates supercritically from 
the layer state. Therefore, close to the bifurcation point, 
this weakly phase separated steady state can be found by 
linearizing Eq. ([5]) around the layer state u — U (x) . In 
the time-dependent hydrodynamic framework, this lin- 
ear analysis corresponds to a marginal stability analysis 
of the layer state with respect to small perturbations in 
the y- and z-directions. 

Substituting u{x,y,z) = U(x) + ipk(x) cos k y y cos k z z 
into Eq. ([5]) and linearizing with respect to the small 
correction ipk(x) cos k y y cos k z z, we obtain the following 



linear equation for the new function 4>k(x) = F ipk{x): 

^(z)-(fc 2 + ^)^(x)=0. (9) 

Here k 2 = k 2 + k 2 , the functions F and Q are evaluated 
at u = U{x), and Qu denotes the u derivative of Q(u) 
evaluated at u = U(x). The boundary conditions for 
Eq. © are 

<j)' k {x = 0) = and 4> k {x = 1) = . (10) 

Equation can be interpreted as a Schrodinger equa- 
tion for a single particle in a one-dimensional potential 
V(x) — AQu/F, while the quantity — k 2 serves as the 
(negative or zero) energy eigenvalue. We solved this 
eigenvalue problem numerically for different A and /. 
Figure [4] shows the resulting marginal stability curves 
£;* = &*(/) for four different values of A > A c . Assuming 
that the instability is non-oscillatory at the onset, one can 
interpret the marginal stability results as linear instabil- 
ity of the layer state below the corresponding curve, and 
linear stability above the curve. The instability is possi- 
ble only within the spinodal balloon: the borders fi and 
fi of the spinodal interval correspond, at fixed A, to zero 
eigenvalues k — > 0. 
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FIG. 4: Marginal stability curves for four different values of 
A are plotted in rescaled coordinates: fcA _1//2 versus /A 1 / 2 . 
For a fixed A the layer state is linearly stable above the cor- 
responding curve and unstable below the curve. 

It can be seen in Fig. [5] that, in the rescaled coordi- 
nates, the marginal stability curves for different A orig- 
inate (almost) at the same point of the horizontal axis 
/A 1 / 2 . Furthermore, the maxima of all the curves are 
almost equal. Like in the 2D case, the first property can 
be explained analytically by considering the dilute limit 
of the problem, while the second property results from 
the strong localization of the eigenfunctions <pk(x) near 
the elastic wall [12} . 

Having found the eigenvalues fc* (/, A) , we can deter- 
mine the critical (minimum) lateral aspect ratios for a 
phase separation. Indeed, the zero heat flux conditions 
at the walls y — 0, y — A y , z — and z — A z (recall that 
we are using rescaled coordinates) yield the quantization 
rules k y = (n y Tr)/A y and k z = (n z Tr)/ A z , where n y ,n z = 



5 



0, 1,2,.. .. Therefore, the minimum value of A y (A z ) for 
a phase separation only in the y-direction (correspond- 
ingly, only in the z-direction) is A y = A c z = 7r/fc*(/, A). 
For example, for A = 8 x 10 3 and / = 0.011 we obtain 
kl ~ 21.0, therefore A y = A% = n/K ~ 0.69. In order to 
have a phase separation in both lateral directions y and 
z, the aspect ratios A y and A z must obey the inequality 

J_ 1 fcg(/,A) 

III. MD SIMULATIONS 
A. Method 

We performed a series of event-driven MD simulations 
of this 3D system using a standard algorithm described 
by Rapaport [24|. Simulations involved N hard spheres 
of diameter d — 1 and mass m = 1. After each collision 
of particle i with particle j, their relative velocity was 
updated according to 

v-j =Vij-(l + r) (vij ■ r,,)r,, , (11) 

where is the precollisional relative velocity and f-y = 
rij/\fij\ is a unit vector connecting the centers of the 
two particles. The "thermal" wall was kept at constant 
temperature To that we set to unity. We used a stan- 
dard thermal wall implementation, see e.g. Ref. [25| . 
p. 173-177. Particle collisions with the rest of the walls 
were assumed elastic. The natural time unit of the MD 
simulations is dim/To) 1 ^ 2 — 1. The initial spatial distri- 
bution of (non-overlapping) particles was uniform in all 
simulations. The initial particle velocity distribution was 
uniform in the direction angles, while the absolute value 
vq of the velocity of each particle was chosen to be such 
that (1/2) mv 2 = (3/2)(l - l/iV)T . In all simulations 
the velocity of the center of mass of the particles at t = 
was zero. That the transients died out and the system 
reached a steady state was monitored by (i) measuring 
the total energy of all particles versus time, and (ii) mea- 
suring the coordinates of the center of mass versus time. 

As a test simulation, we performed a simulation of 
a system which dimensionless parameters A and / are 
within the spinodal balloon but which still cannot phase 
separate because of too small lateral dimensions, see 
Fig. [TJ The inset shows a late-time snapshot of the 
system as observed in the MD simulation. As one can 
see from Fig. [TJ the measured particle number density, 
rescaled to function of the rescaled distance from 

the driving wall is in good agreement with our hydro- 
static calculations. 



B. Phase separation 

Remaining within the spinodal balloon, and increasing 
one of the lateral dimensions of the system, we observed 



phase separation as expected from the theory, see Fig. [5] 
A dense cluster develops in one of the two corners of 
this quasi-2D Hele-Shaw cell, at the wall opposite to the 
driving wall. Quantitative diagnostics are provided by 
the plots of the three center-of-mass coordinates of the 
system versus time, shown in the lower panel of Fig. [5] 
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FIG. 5: Upper panel: a late-time snapshot of an MD simula- 
tion with N = 10 s particles of r = 0.89945 in the container 
with dimensions L x = L y = 500d and L z — 50d. The upper 
wall is the driving wall. The dimensionless parameters of the 
system A = 5 x 10 4 and / = 0.0057 correspond to a point 
within the spinodal balloon of Fig. [2] Lower panel: the ab- 
solute values of the three center-of-mass (COM) coordinates 
(measured in the units of the particle diameter d) versus time 
[measured in the units of d^m/To) 1 ^ 2 ]. The center of the con- 
tainer is at the origin. 



Figure [6] shows another example of phase separation 
for A and / within the spinodal balloon, but this time in 
the case when both lateral dimensions L y and L z are suf- 
ficiently large. As one can see, a dense stripe-like cluster 
forms along one of the edges of the wall opposite to the 
driving wall. 

Both Fig. [5l and Fig. [B] show phase separated states 
with a 2D, rather than 3D, structure. An example of a 
fully 3D structure is shown in Fig. [TJ Here a fully 3D 
dense cluster (a "drop" ) develops in one of the corners at 
the wall opposite to the driving wall. 
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FIG. 6: Upper panel: a late-time snapshot of an MD sim- 
ulation with N = 5 x 10 5 particles of r = 0.899452 in the 
container with dimensions L x — 200d and L y — L z = 400d. 
The upper wall is the driving wall. The dimensionless param- 
eters of the system A = 8 x 10 3 and / = 0.011 correspond to 
a point within the spinodal balloon of Fig. [2] Lower panel: 
same as in Fig. 
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FIG. 7: Upper panel: a late-time snapshot of an MD sim- 
ulation with N = 340736 particles of r = 0.899452 in the 



container with dimensions L x 



200d and L v 



400d. 



The upper wall is the driving wall. The dimensionless pa- 
rameters of the system A = 8 x 10 s and / = 7.529 x 10~ 3 
correspond to a point within the spinodal balloon of Fig. [2] 
Lower panel: same as in Fig. [5] 



IV. DISCUSSION 

As we have shown, granular hydrodynamics predicts 
negative lateral compressibility and, therefore, phase- 
separation instability of the basic state of a granulate 
driven by a thermal wall of a rectangular container at 
zero gravity. When the lateral dimensions of the con- 
tainer are sufficiently large, the negative compressibility 
causes a van der Waals-like phase separation instability. 

Our MD simulations are in agreement with hydrostatic 
predictions (of course, if we disregard small fluctuations 
caused by the discreteness of the particles). In the lan- 
guage of hydrostatics, a broken-symmetry steady state is 
described by either a 2D, or a fully 3D solution of the 
nonlinear partial differential equation (|5|) subject to the 
fixed mass constraint J7]) and the boundary conditions. 
Such steady-state solutions can be obtained only numer- 
ically (see Ref. @ for 2D examples). Because of the 
translational symmetry of the steady-state equations in 
the y and z directions, bounded solutions satisfying the 
no- flux boundary conditions in these directions, must be 
either independent of the y and z coordinates, or pe- 
riodic in them. Furthermore, by analogy with the 2D 
setting, one should expect that dynamic coarsening se- 



lects a periodic steady state solution with a maximum 
spatial period (equal to twice the container size in the 
corresponding direction [Hj]). 

When one of the lateral aspect ratios, say A y , is larger 
than the critical value A C (A, /), while the other one, A z , 
is smaller than A c , a 2D pattern should develop, and 
this is indeed what we observed in Fig. When both of 
the lateral dimensions are larger than A c , we observed 
dense clusters of two different morphologies: either a 2D 
morphology, like the one shown in Fig. [6l or a fully 3D 
morphology, like the one shown in Fig. [7J 

What will happen when the parameters A and / are 
within the spinodal or binodal balloons, and its lateral 
dimensions L x and L y are much larger than Ac? We 
expect that multiple "drops" (or bubbles) will nucleate 
at the wall opposite to the driving wall and undergo dy- 
namic coarsening, qualitatively similar to Ostwald ripen- 
ing [2(| , before reaching the final state with a single drop 
(or bubble). The Ostwald ripening regime is beyond the 
reach of our present computing resources. 

When the lateral aspect ratios A y and A z are just 
above (or just below) the critical value A c = A C (A, /), 
the phase separation, as predicted by hydrostatics, 
should be "weak" and look, in the supercritical case, 
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as a small-amplitude modulation of the layer state 
d EL HI HI Q HI. In 2D such a system experiences 
large fluctuations [l4j , and it would be interesting to find 
out whether large fluctuations persist in 3D. 

We also performed a series of MD simulations for more 
realistic conditions, using a (truly) vibrating wall instead 
of the thermal wall, and allowing for inelastic particle 
collisions with the walls. Qualitatively, the results have 
not changed: for sufficiently high inelasticity of particle- 
particle collisions we observed phase separation for in- 
termediate values of the volume fraction, and no phase 
separation for too a small or too a large volume fraction. 

Our theory and simulations assumed a a zero gravity. 
To what extent is the van der Waals-like phase transition 
sensitive to the presence of a small gravity, with acceler- 
ation g, directed towards the thermal wall? In a 2D set- 
ting this question was addressed by Khain and Meerson 
|27j . Assuming dilute limit of granular hydrodynamics, 
they found that, as the Froude number F = mgL x /T 
increases, the phase separation crosses over to "thermal" 
granular convection. One can expect a similar scenario 



in 3D as well, though this question has not been yet been 
considered in detail. 

In summary, the van der Waals-like phase separation in 
2D and 3D provides a useful and rich prototypical model 
system for testing the ideas and methods of granular dy- 
namics. By focusing our attention on 3D in this work, 
we broke ground for a future investigation of this fasci- 
nating phase separation phenomenon in reduced gravity 
experiments. 
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